Microfluidic capturing-dynamics of paramagnetic bead suspensions 
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We study theoretically the capturing of paramagnetic beads by a magnetic field gradient in a 
microfluidic channel treating the beads as a continuum. Bead motion is affected by both fluidic and 
magnetic forces. The transfer of momentum from beads to the fluid creates an effective bead-bead 
interaction that greatly aids capturing. We demonstrate that for a given inlet flow speed a critical 
density of beads exists above which complete capturing takes place. 

PACS numbers: 47.15.Pn, 47.55.Kf, 47.60.+i, 41.20.-q 
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I. INTRODUCTION 

Recently, there has been an increasing interest in using 
magnetic beads in separation of, say, biochemical species 
in microfluidic systems P, Q] . The principle is to have 
biochemically functionalizcd polymer beads with inclu- 
sions of superparamagnetic nanometersize particles of, 
for example, magnetite or maghemite. They attach to 
particular biochemical species and can be separated out 
from solution by applying external magnetic fields. As 
most biological material is either diamagnetic or weakly 
paramagnetic, this separation is specific. Paramagnetic 
particles in fluids are also used to measure the suscepti- 
bility of, for example, magnetically labelled cells by mea- 
suring particle capture or motion in a known field 0, Q . 

In this paper we study microfluidic capturing of para- 
magnetic beads from suspension by modeling the beads 
as a continuous distribution [jj. The separation of sus- 
pended paramagnetic beads from their host fluid is an 
important process as it decides operating characteristics 
for practical microfluidic devices. It involves an interplay 
between forces of several kinds governing the dynamics of 
the process: (a) Magnetic forces from the application of 
strong magnetic fields and field gradients, (b) Drag forces 
due to the motion of the beads with respect to the host 
fluid, (c) The trivial effect of gravity, which we ignore in 
the following. We emphasize the effects of bead motion 
on the fluid flow as this gives rise to a hydrodynamic 
interaction between the beads. As we have noted in a 
previous few-bead study, this interaction is more impor- 
tant than the magnetic bead-bead interactions 6]. It is 
created by drag forces in two steps: First, drag transfers 
momentum to the fluid from the beads moving under the 
influence of external forces. Second, the modified flow 
changes the drag on and thus motion of other beads. 



II. MODEL 

As sketched in Fig. [Ha) we consider a viscous fluid 
(water) flowing in the x direction between a pair of par- 
allel, infinite, planar walls. The walls are placed parallel 
to the xy plane at z — and z — h, respectively. A 
steep magnetic field gradient is generated by a parallel 



pair of closely spaced, infinitely long, thin wires along the 
y direction separated by d and carrying opposite currents 
±1. The system is translation invariant in the y direction 
thereby reducing the simulation to a tractable problem in 
2D. The simulation domain is defined by < x < L and 
< z < h with L = 350 /im and h = 50 /im. The wires 
intersect the xz plane near (x, z) = (250 /im, 55 /im) just 
above the top plate. Paramagnetic beads in suspension 
are injected into the microfluidic channel by the fluid flow 
at x — 0. They are either exiting the channel at x = L 
or getting collected at the channel wall near the wires. 

When a suspension of beads is viewed on a sufficiently 
large scale compared to the single bead radius a but on 
a scale comparable to density variations, we can describe 
the distribution of beads in terms of a continuous, spa- 
tially varying bead number density c. We consider a sus- 
pension of beads with radius a — 1 /im and denote the 
initial number density at x = by c . The four basic 
constituents of the model are described in the following. 
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r\ = 1 mPa s 
q = 10 3 kg/m 3 
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Id = 2 /iAm 
c = 10 13 m~ 3 



Id = 12 /iAm 
c = 10 13 m~ 3 



Id = 40 /iAm 
c = 10 13 m~ 3 




FIG. 1: (a) Sketch of the microfluidic system with L = 350 fim 
and h = 50 /im. A suspension of paramagnetic beads enters at 
x — with a parabolic Poiseuille flow profile, u , and leaves at 
x = L. Beads are caught by the pair of wires placed 100 /im 
from the outlet at the top and carrying currents ±1. (b)- 
(d) Simulated stationary density of the beads ranging from 
zero (white) to c (black) for increasing values of the current- 
distance product Id as indicated. At x = the concentration 
is c = 10 13 m~ 3 and the maximum flow speed is 300 /jm/s. 
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Magnetic force. The beads are paramagnetic with a 
magnetic susceptibility % = 1. In an external magnetic 
field H cxt (r) the force on such a bead is 

" cxt — f^o 

/ (M • V)H cxt dV 

J bead 

= 4wi 3^_(H ext - V)H ext (1) 

assuming that the bead is so small that we can take the 
external field H oxt to be approximately constant across 
the bead, i.e., a|VH cxt | <C |H cxt | when determining the 
magnetization M. 

As mentioned, H cxt in this study arises from a pair 
of current carrying wires. It is determined in the fol- 
lowing manner. From Ampere's law, we readily find 
the magnetic field, H, around a straight circular wire, 
H(r) = J x r/(27rr 2 ), where the electrical current vec- 
tor J is along the wire orthogonal to the position vector 
r which is in the xz plane. The magnetic field from the 
two closely spaced anti-parallel wires is found by decreas- 
ing the separation d and increasing the current, I = |J|, 
while keeping the product Id constant, 

This together with Eq. yields 

F C xt = --^/ioa 3 (/c0 2 4; ( 3 ) 
7r x + o r° 

a manifestly attractive central force (from the mid-point 
of the wires), independent of the direction of d. 

Fluid flow. The beads are suspended in a fluid of vis- 
cosity r] and density g that is launched at x = with a 
parabolic velocity profile, u Q , and flows past the wires. In 
microfluidics inertial effects are unimportant compared 
to drag, so the small beads in suspension almost always 
move with constant velocity relative to the fluid. Ex- 
cept for acceleration phases shorter than microseconds 
the external forces are exactly balanced by drag Q ■ The 
momentum transfer from beads to fluid is included by 
adding a bulk force term, cF ext (r) to the Navier-Stokes 
equation. This bulk force term is proportional to the 
number density c of beads and the magnetic external 
force F oxt on an individual bead at position r. The ve- 
locity u of the fluid is given by 

gd t u + q(u- V)u = -Vp + 77V 2 u + cF cxt , (4) 

along with the incompressibility condition V • u = 0. 

Bead motion. To complete the set of equations, it is 
necessary to have an equation of motion for the bead 
number density c. As the beads neither appear nor dis- 
appear in the bulk, c must obey a continuity equation 

9 t c + V-j = 0, (5) 




FIG. 2: The fraction f3 of beads caught as function of the 
current-distance product Id for twenty different flow speeds 
(50 - 1000 /im/s; indicated by the arrows). Larger current 
leads to higher /3; faster flow to smaller /3. In this simulation 
the initial concentration is low, c = 10 13 m -3 . Inset: Rate 
7 ca p of bead capture as function of Id, for the flows above. 
The faster flow or the larger current, the higher 7 cap . 



where the bead current j is defined by the Nernst-Planck 
equation Q 

j = -DVc + cu + c6F 0Xt (6) 

with diffusivity D and bead mobility b — l/(6wqa). 

For our spherical beads the diffusivity is given by the 
Einstein expression D — kT / (Qirrja) which for water at 
room temperature equals 2.2 x 1CP 13 m 2 /s. In the simu- 
lations below, however, we artificially increase the mag- 
nitude of D in order to stabilize the computations and to 
use a coarser mesh and thus save computation time. 

Boundary conditions. In addition to the bulk equa- 
tions , ||SJ|, and (JSJ), we need appropriate boundary 
conditions. As the beads move out to the walls of the 
domain and settle there, merely demanding that the nor- 
mal component of the bead current vector j vanishes is 
not correct, rather, it must be free to take on any value 
as long as it is directed into the wall. As beads do not en- 
ter the bulk from the walls (by assumption once settled, 
beads stick) we demand that the normal current compo- 
nent is never directed into the liquid. For the fluid we 
demand the usual no-slip condition at the walls. 

At the inlet x = of the microfluidic channel we as- 
sume that the fluid comes in with the constant initial 
number density c and with a parabolic fluid velocity 
profile u Q . At the outlet x = L we let the bead current 
take on any value, while the fluid pressure is zero. 



3 




4 [A 2 /xms] 



FIG. 3: The fraction /3 of beads caught versus (Id) 2 /u , the 
ratio of the current-distance product squared and the fluid 
flow velocity. This demonstrates scaling in the competition 
between capturing and flushing: the twenty curves from Fig. H 
approximately collapse to one single. 
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FIG. 4: Fraction (3 of beads caught as function of initial 
bead number density with and without the bulk force term 
cF oxt in Eq. (gj. The fixed values for the current-distance 
product Id and the maximum in-flow speed u are shown. At 
low densities less than 50% are caught; at high densities the 
collective motion of the beads leads to 100% capture. 



III. RESULTS 

Having set up the equations for bead and fluid motion, 
they are solved with the finite element method on a mesh 
with ~ 10 4 elements refined in the vicinity of the wires. 
To this end we employ the finite element solver software 
package Femlab 9]. The parameter values for the fluid 
are those of water, rj = 1 mPas and g — 10 3 kg/m 3 , while 
for the beads a = 1 fim and c Q = 10 13 to 10 16 m -3 . 

To study capturing we must keep track of which beads 
are captured and which are flushed through the channel 
with the flow. This is done by calculating the rates 7* 
by which the beads are either captured or transported 
in/out at each of the four boundary segments i of the 
channel (inlet, outlet, upper wall, and lower wall). By 
integration of the normal components of the bead cur- 
rents along each segment i, we find 

ji= fi-ndii. (7) 

J i 

The rate of capture is 7 cap = 7i OWO r + 7u PP cr- In steady 
state the conservation of beads enforces 7i n iot + 7cap + 
7outiot = 0, which provides a useful check of the simu- 
lation results. The primary control parameters are the 
current-wire distance product Id, the maximum fluid in- 
flow speed u a , and the bead number density, c Q . The 
product Id decides the magnetic force which captures the 
beads against the fluid flow. As we are investigating ef- 
fects of bead-bead interaction, our interest is properties 
that depend on the bead number density, in particular 
those that do so nonlinearly. 

Electrical current and fluid flow. The effects of having 
electrical wires near, and thus a magnetic field gradient 
in, the channel is illustrated in Fig. IHb)-(d) for three 
values of the current-distance product Id. At small val- 
ues of Id only a narrow region is emptied but increasing 



the current the region expands until it covers the width 
of the channel. 

A simple measure of the capturing is the ratio (3 of the 
bead capture rate 7 cap to the bead in-flow rate 7i n i e tj 

^ _ "capture rate" _ 7 cap , g s 
"in-flow rate" 7iniet 

If capturing dominates (3 tends to unity, if flushing dom- 
inates (3 tends to zero. Fig. [21 shows this in that slow 
flow and strong current leads to a high j3 whereas fast 
flow and weak current leads to a small value. The rate 
7cap of bead capture as function of wire current and flow 
velocity is illustrated in the inset of Fig. [21 

If there is a competition between magnetic capturing 
and flushing, then we expect that the data can be de- 
scribed essentially by the ratio of the magnetic forces to 
the inlet fluid flow speed u a . The force is proportional to 
the square of the current-distance product Id. In Fig. [31 
we plot the data from Fig. [21 as function of (Id) 2 /u a and 
see that the data mostly collapses onto a single curve. 
The collapse is not perfect and is not expected to be as 
the underlying flow and bead distribution patterns (see 
Fig.[TJ arc different for different flows and magnetic fields. 

Interactions and concentration. The second point we 
wish to make is that modification of the overall flow, and 
the effective bead-bead interaction this entails, is signifi- 
cant for bead capturing. We can study the effect by ex- 
cluding momentum transfer to the fluid flow due to the 
bulk force term cF 0Xt in the Navier-Stokes equation 
At high bead number densities the force acting on the 
beads contributes a significant force on the fluid affect- 
ing fluid flow and spawning the effective interaction. The 
strength of this interaction must thus depend on the den- 
sity of particles. This is illustrated in Fig. capturing 
was simulated at fixed in-flow speed, u Q — 300 /im/s, 
and a fixed value of the current- wire distance product, 
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A/3 = /3i nc l — /3 ex cl 








FIG. 5: The difference A/3 = /3mci — /3 C xci in captured bead 
fractions between two situations: including and excluding the 
bulk force term cF oxt in the Navier-Stokes equation (QJ. At 
high concentrations (c > 10 15 m -3 ) there is an apprecia- 
ble difference between including and excluding the bulk force 
term, corresponding to hydrodynamic bead-bead interactions. 

Id = 8 /iAm, but for varying bead number densities c 
ranging from 10 13 to 10 16 m -3 . At low densities we find 
that capturing is roughly independent of density and the 
fraction (3 of beads captured has some intermediate value, 
however, for high densities all beads are caught. In con- 
trast, leaving out the bulk force term cF cxt in the Navier- 
Stokes equation, i.e., the force acting on the fluid, gives 
concentration independence as shown in Fig. 0] 

As can be seen in Fig. [5] a complementary way of ex- 
hibiting the importance of the bulk force term is to plot 
the difference A(3 — /3 inc i — /3 oxcl between including and 
excluding cF cxt as a function of concentration and the 
current-distance product. This shows that interactions 
makes an appreciable difference at high concentrations 
and intermediate magnetic fields. 

Diffusion constant. Even for the small beads of radius 
a = 1 /im, the diffusion constant given by the Einstein 
relation is small compared to the dimensions entering the 
problem. The time-scale for a bead to diffuse across the 
channel is Tdiff ~ h 2 /D. If we are to see the influence of 
diffusion competing with bead advection, then the rele- 
vant quantity is the Peclet number hu/D which is advec- 
tion time-scale r a d v ~ hju over the diffusion time-scale. 



When this number is larger than unity, which it is except 
for artificially large diffusion constants, then convection 
dominates. In the simulations the diffusion constant is 
increased artificially up to 10~ n m 2 /s in order to help 
numerical convergence. But we have verified that values 
smaller than 10 -10 m 2 /s have no influence on the results. 



IV. DISCUSSION AND CONCLUSION 

We have studied microfluidic capture of paramagnetic 
beads in suspension. The three main findings of work are: 
the approximate scaling shown in Fig. El the existence of 
a critical bead density for capture shown in Fig. and 
the qualitative difference for capturing between models 
with and without the hydrodynamic bead-bead interac- 
tion shown in Figs. |3] and 

Clearly, it is very important for the capture process to 
include the action of the beads on the host fluid medium. 
Simply leaving it out can give qualitatively wrong results 
for high concentrations of beads. This casts some doubt 
on the measurement of cell susceptibility through captur- 
ing as it depends on cell concentration 0, 0] . Deduction 
of susceptibilities from single bead or cell considerations 
together with measurements at high bead or cell concen- 
tration is suspect. Care must be taken to compare with 
standards of known and similar susceptibility, size, and 
concentration. 

The effective bead-bead interaction greatly helps cap- 
turing. It should make detectable differences depending 
on whether there are a few or hundreds of particles in a 
channel at a time in actual experiments especially when 
the flow and magnetic field are such that the beads are 
barely caught one by one. This interaction should be 
considered when choosing operating conditions for mi- 
crofluidic devices based on capturing of beads as higher 
bead number densities potentially eases requirements for 
external magnets and allows faster flushing. We hope 
that experimental studies will be initiated to verify this 
prediction of our work. 
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